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Preface 


The  study  documented  herein  was  imdertaken  as  part  of  the  Work  Unit  “Innovative  Geophysical 
Technologies  for  Enhanced  Buried  Unexploded  Ordnance  (UXO)  Discrimination”  (AF25,  6.2)  funded  by 
the  U.S.  Army  Engineer  Research  and  Development  Center  (ERDC),  under  the  Enviromnental  Quality 
Technology  (EQT)  Program.  Dr.  Dwain  K.  Butler,  Geotechnical  and  Structures  Laboratory  (GSL),  is 
Principal  Investigator  for  the  work  unit,  and  Dr.  John  M.  Cullinane,  Environmental  Laboratory,  is  EQT 
Program  Manager.  The  work  was  performed  at  the  University  of  British  Columbia  (UBC),  Vancouver, 
British  Columbia,  Canada,  and  at  the  ERDC,  Vicksburg,  MS,  facility.  Funding  of  the  UBC  effort  was  by 
means  of  Army  Research  Office  Grant  41262-EV,  “Locating  and  Characterizing  Unex^ploded  Ordnance 
Using  Time  Domain  Electromagnetic  Induction.”  Dr.  Russell  Harmon,  Terrestrial  Sciences  Program 
Manager,  administered  the  UBC  grant. 

The  work  was  performed  dining  the  period  1  April  2000  and  3 1  March  2001  by  Mr.  Leonard  R.  Pasion 
and  Dr.  Douglas  W.  Oldenburg,  UBC,  with  the  assistance  of  Mr.  Jose  L.  Llopis  and  Dr.  Butler  for  the  work 
at  the  ERDC  facility.  The  work  is  based  on  a  technique  for  interpreting  time  domain  electromagnetic 
(TEM)  induction  data  proposed  by  Pasion  (1999),  as  listed  in  the  Reference  section.  Part  I  of  the  report 
documents  development  of  a  discrimination  algorithm  for  UXO  based  on  TEM  induction  and  validates  the 
algorithm  by  application  to  a  Geonics  EM63  TEM  data  set  obtained  at  the  ERDC  UXO  Test  Site, 
Vicksburg,  MS  (Pasion  and  Oldenburg  2001).  Part  II  of  the  report  contains  descriptions  of  the  extensive 
Geonics  EM63  TEM  signature  database  obtained  at  the  ERDC  UXO  test  site  during  the  study.  This  part 
outlines  the  tests  verifying  the  applicability  of  the  two-dipole  model  and  concludes  with  the  discrimination 
algorithm  being  applied  to  a  number  of  UXO  targets. 

At  the  time  of  publication  of  this  report.  Dr.  Michael  J.  O’Connor  was  Director,  GSL,  and  Dr.  James  R. 
Houston  was  Director,  ERDC.  Commander  and  Executive  Director  was  COL  John  W.  Morris  HI,  EN. 
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Parti 


A  Discrimination  Algorithm  for  Unexploded 
Ordnance  (UXO)  Using  Time  Domain 
Electromagnetics  (TEM) 


1  Introduction 


An  explosive  ordnance  is  a  munition  that  is  either  launched  or  fired  with  the  intent  of  detonation  at  a  specified 
target.  An  unexploded  ordnance  (UXO)  is  an  explosive  ordnance  that,  because  of  some  malfunction,  remains 
undetonated.  As  a  result,  the  ordnance  can  be  found  at  the  ground  surface,  partially  buried,  or  buried  at  a  depth 
of  up  to  8  m  beneath  the  surface.  The  remediation  of  UXO-contaminated  land  has  been  made  a  high  priority  by 
the  United  States  Department  of  Defense  in  order  to  either  maintain  safe  usage  for  continuing  military  operations 
or  to  permit  land  transfer  to  the  private  sector.  Practical  and  cost-effective  strategies  for  remediation  require  both 
detection  of  possible  targets  and  the  ability  to  discriminate  between  UXO  and  contaminating  scrap  metal. 

The  detection  of  buried  metallic  objects  can  be  accomplished  with  a  variety  of  geophysical  sensing  techniques. 
Time  domain  electromagnetic  (TEM)  induction  surveys  have  been  successful  in  detecting  both  ferrous  and  non- 
ferrous  metaUic  objects  near  the  soil  surface  and  are  a  mainstay  among  technologies  currently  utilized  in  UXO 
clearance  projects.  In  the  TEM  method,  a  time  varying  magnetic  field  is  used  to  illuminate  a  conducting  target. 
This  primary  field  induces  surface  currents  on  the  target  which  then  generate  a  secondary  magnetic  field  that  can 
be  sensed  above  ground.  With  time,  the  surface  currents  diffuse  inward,  and  the  observed  secondary  field  conse¬ 
quently  decays.  The  rate  of  decay,  and  the  spatial  behavior  of  the  secondary  field,  are  determined  by  the  target’s 
conductivity,  magnetic  permeabihty,  shape,  and  size. 


Identification  of  a  UXO  from  electromagnetic  sensor  data  remains  a  major  hurdle  in  reducing  the  high  costs  of 
remediation  projects.  It  has  been  reported  that  approximately  70  percent  of  remediation  costs  are  currently  being 
used  to  excavate  nonordnance  items  (Butler  et  al.  1998).  The  development  of  discrimination  algorithms  can  be 
roughly  categorized  as  either  model-based  or  data-based.  Data-based  algorithms  are  pattern  recognition  proce¬ 
dures  that  compare  a  library  of  catalogued  responses  from  various  UXO  items  to  measured  responses  (Damarla 
and  Ressler  2(KX)).  Model-based  algorithms  use  either  an  exact  or  approximate  forward  modeling  algorithm  to  de¬ 
termine  a  set  of  model  parameters  needed  to  replicate  the  measured  responses,  and  subsequently  relating  the  model 
parameters  to  physical  parameters  (Khadr,  Barrow,  and  Bell  1998).  One  such  model-based  techmque  that  has  been 
the  focus  of  much  recent  research  is  the  determination  of  the  time  constants  of  the  TEM  response,  or  equivalently, 
the  poles  of  the  frequency  domain  signal,  to  identify  the  buried  target  (Snyder  et  al.  1999;  Baum  1997 ;  Collins  et 
al.  1999).  A  method  that  represents  a  hybrid  of  the  model-based  and  data-based  algorithm  is  under  development  at 
Blackhawk  Geometries  (Grimm  2()(X)).  In  that  approach,  a  spheroid  modeller,  working  jointly  with  a  model-based 
inversion  algorithm,  generates  a  library  of  model  parameters  which  can  then  be  operated  upon  by  a  neural  network 
classifier  for  comparison  with  parameters  derived  from  the  raw  sensor  signal. 


In  this  paper  we  present  a  model-based  TEM  data  interpretation  algorithm  which  estimates  the  basic  shape  (rod¬ 
like  or  plate-like)  and  magnetic  character  (ferrous  or  nonferrous)  of  a  buried  metallic  object.  We  first  present  an 
approximate,  forward  model  that  represents  the  time  domain  response  of  a  metallic  object  as  a  pair  of  perpendicular 
dipoles  located  at  the  center  of  the  buried  target.  This  form  of  model  was  suggested  to  us  in  a  personal  communi¬ 
cation  from  J.D.  McNeill  (1997).  The  strengths  of  these  dipoles  decay  with  time,  and  the  parameters  that  govern 
the  time  decay  behavior  are  related  to  the  conductivity,  permeability,  shape,  and  size  of  the  buried  target.  The 
parameterization  is  simple  and,  thus,  convenient  to  use  in  data  fitting  procedures.  We  next  describe  an  inversion 
scheme  to  recover  the  model  parameters  from  TEM  data.  Since  these  parameters  encapsulate  information  about 
the  physical  attributes  of  the  target,  we  can  attempt  to  use  them  to  determine  if  the  target  is  ferrous  and  if  the 
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geometry  is  rod-like  (most  likely  a  UXO)  or  plate-like  (most  likely  a  nonordnance  item).  Empirical  relationships 
are  developed  that  link  the  model  parameters  to  the  physical  parameters  of  the  target  and  these  relationships  form 
the  basis  of  our  algorithm.  We  conclude  with  the  application  of  the  algorithm  to  a  synthetic  data  set  contaminated 
with  noise,  along  with  field  data  sets  taken  over  a  buried  UXO  and  buried  metallic  scrap. 
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Chapter  1 


Introduction 


2  Development  of  Approximate  Forward 
Modelling 


In  order  to  invert  measured  TEM  data  for  the  physical  parameters  of  the  target,  it  is  necessary  to  have  a  forward 
model  to  describe  the  TEM  response  for  a  buried  metallic  object.  We  can  restrict  our  search  for  response  solutions 
to  axi-symmetric  metallic  targets,  since  this  geometric  subset  adequately  describes  aU  UXO  and  the  majority  of 
buried  metallic  scrap  encountered  in  a  remediation  surv^.  Unfortunately,  analytic  expressions  for  the  time  domain 
response  are  restricted  to  a  metallic  sphere,  and  even  an  expression  for  a  permeable  and  conducting  nonspherical 
axi-symmetric  body  is  not  available.  Numerical  solutions  of  Maxwell’s  equations,  under  continual  development, 
are  promising  (e.g.  Haber,  Ascher,  and  Oldenburg  2000;  Carin  2000;  Hiptmair  1998);  however,  the  computational 
time  requirements  for  obtaining  a  solution  still  make  them  impractical  for  use  as  part  of  a  rigorous  inversion 
procedure.  Our  approach,  therefore,  is  to  use  an  approximate  forward  modelling  that  can  adequately  reproduce 
the  measured  electromagnetic  anomaly  in  a  minimal  amount  of  time.  The  validity  of  this  reduced  modelling  stiU 
needs  further  testing,  but  the  empirical  tests  carried  out  here  suggest  that  it  can  be  useful  in  practice. 

The  development  of  the  approximate  forward  modelling  is  presented  in  four  steps.  We  begin  with  the  response 
of  a  sphere,  so  that  the  magnetic  polarization  dyadic  M  is  introduced.  This  dyadic  is  then  altered  so  that  it  is 
applicable  to  an  axi-symmetric  body.  This  generates  the  ”two-dipole”  model  mathematically.  Next,  we  introduce 
a  parameterization  for  the  time  decays  of  each  of  the  two  dipoles,  and  finally,  we  combine  everything  to  generate 
our  approximate  forward  modelling. 


Response  for  a  Spherical  Body 

Consider  a  permeable  and  conducting  sphere  of  radius  a  illuminated  by  a  uniform  primary  field  B  ^  (Figure  la). 
At  a  time  t  =  0  the  primary  field  is  terminated,  and  eddy  currents  are  induced  in  the  sphere;  they  subsequently 
decay  because  of  the  finite  conductivity  of  the  sphere.  The  secondary  field  B  generated  by  the  decaying  currents 
is  dipolar: 

(1) 

where  m  (t)  is  the  dipole  moment  induced  at  the  center  of  the  sphere  at  time  t,  r  is  the  distance  between  the 
observation  point  and  the  sphere  center,  f  is  the  unit  vector  pointing  from  the  sphere  center  to  the  observation  point 
P,  and  i  is  the  identity  dyadic.  The  dipole  moment  is 

m  (t)  =  —  bP  T®  (t)  (2) 

1^0 


where 


CXD 

(t)  = 

s=l 


—  Qsi 

e  r 


Qs  it^r  1)  (/^r  2) 


(3) 
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where  r  =  afia^,  and  Hr  =  m/mo  is  the  relative  permeability  (Kaufmann  and  Keller  1985).  In  general,  the 
magnetic  permeability  of  highly  permeable  materials  is  a  function  of  many  parameters,  including  the  strength  of 
the  incident  magnetic  field,  temperature,  and  magnetic  history.  However,  calculated  TEM  responses  assuming 
a  constant  permeability  of  =  150  for  steel  and  =  1  for  aluminum  compared  well  with  laboratory  TEM 
measurements  of  steel  and  aluminum  targets  (Pasion  1999).  Therefore  we  feel  that  Equations  1  to  3  are  suitable 
for  the  analysis  that  follows  here.  The  values  Qs  roots  to  the  transcendental  equation 


tangfg  = 


(/^r  1)  Qs 

ql  +  [Hr  -  1)  ■ 


(4) 


Equations  1  to  4  reveal  that  the  secondary  B-field  of  a  sphere  in  a  uniform  primary  field  is  equivalent  to  the 
B-field  of  a  single  magnetic  dipole  located  at  the  center  of  the  sphere  and  oriented  parallel  to  the  primary  field.  For 
convenience,  we  write  the  relationship  between  the  induced  dipole  and  the  primary  field  as  m  =  M  •  B^,  where 
M  is  the  magnetic  polarizability  dyadic.  For  a  sphere. 


-  27r  -  27r  ^  ^ 

M  =  — L®  (f)  I  =  —  0  (t)  0 

[  0  0  L®(f) 


(5) 


Baum  (1999)  details  the  characteristics  of  the  magnetic  polarizability  dyadic  and  notes  that  the  triple  degeneracy 
of  the  magnetic  polarizability  dyadic  reflects  the  symmetry  of  the  sphere. 


The  sphere  solution  possesses  several  characteristics  that  we  retain  in  the  formulation  of  our  approximate  solution 
for  an  axi-symmetric  target.  Firstly,  the  secondary  field  resulting  from  the  induced  currents  generated  in  a  sphere, 
illuminated  by  a  uniform,  step-off  primary  field,  is  dipolar  at  all  points  outside  the  sphere.  We  will  also  represent  the 
secondary  field  for  more  general  shapes  as  a  dipolar  field  (Equation  1).  A  dipolar  field  approximation  is  reasonable 
for  any  observation  point  far  enough  away  from  any  localized  current  distribution  (Jackson  1975),  and  it  has  been 
reported  that  for  observation  points  greater  than  1  to  2  times  the  target  length,  a  dipolar  field  assumption  is  adequate 
(Casey  and  Baertlein  1999;  Grimm,  Blohm,  and  Lavely  1997).  Indeed,  higher-order  multipoles  induced  in  a  target 
will  decay  at  early  times  (Grimm,  Blohm,  and  Lavely  1997). 


Secondly,  the  induced  dipole  moment  in  the  center  of  a  sphere  is  given  by  the  dyadic  product  M  ■  B^.  This  form 
indicates  that  the  induced  dipole  is  proportional  to  the  projection  of  the  primary  field  along  the  direction  of  the 
induced  dipole.  The  components  of  M  scale  the  strengths  of  the  dipoles.  The  magnetic  polarizability  dyadic, 
in  the  case  of  the  sphere,  contains  the  function  L  ®  (t)  that  contains  aU  the  information  about  the  time  decay  of 
the  sphere  and  it  depends  upon  the  material  properties,  shape,  and  size  of  the  target.  Our  hypothesis  is  that  more 
general  metallic  shapes  can  also  be  approximately  modell^  with  an  induced  dipole  equal  to  the  dyadic  product 
M  •  B^.  However,  choosing  the  right  functional  form  of  M  will  be  crucial. 


Approximating  the  Magnetic  Polarizability  DyacJic  for  an  Axi-Symmetric 
Body 

Analytic  expressions  for  M  for  the  time  domain  response  of  a  peraeable  and  conducting  nonspherical  axi- 
symmetric  body  are  not  available.  Therefore,  we  base  our  form  of  M  on  the  magnetostatic  polarizability  for  a 
spheroid.  Recall  that,  for  the  time  domain  response  of  a  sphere,  the  structure  of  M  is  identical  to  the  structure 
of  the  magnetostatic  polarizability  dyadic  of  a  sphere.  The  analytic  solution  for  the  magnetostatic  response  of  a 
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magnetic  prolate  spheroid  is  equivalent  to  the  field  of  a  magnetic  dipole  induced  at  the  spheroid  center  (Das  et  al. 
1990): 


^spheroid  mi  +  ma 


=  k,  [(z'  •  BP)  z']  +  fca  [(y'  •  BP)  y'  +  (x'  •  BP)  x']  = 


ifca  0  0 

0  A;2  0 

0  0  fci 


BP 


(6) 


where  ki  and  ko,  are  the  polarizability  constants,  which  are  functions  of  the  conductivity,  permeability,  shape,  and 
size.  Equation  6  reveals  that  the  total  induced  dipole  can  be  written  as  the  sum  of  two  orthogonal  dipoles  m  i  and 
m2.  The  first  dipole  moment  mi  is  parallel  to  the  major  axis  (z'  in  fig.  1(b))  of  the  spheroid,  and  its  strength  is 
proportional  to  the  product  of  the  primary  field  along  that  direction  and  the  polarizability  k  i .  The  second  dipole 
moment  is  perpendicular  to  the  major  axis,  and  its  magnitude  is  proportional  to  the  component  of  the  primary  field 
along  that  direction  and  the  polarizability  k^ .  A  consequence  of  ki  and  k2  being  functions  of  the  spheroid’s  shape 
and  size  is  that  the  orientation  of  the  effective  dipole  will  not  be  solely  determined  by  the  direction  of  the  primary 
field,  as  is  the  case  for  a  sphere.  In  addition,  the  orientation  of  will  be  influenced  by  the  aspect  ratio  of 

the  spheroid. 


The  polarization  dyadic  in  Equation  6  suggests  a  magnetic  polarization  dyadic  for  the  TEM  problem  of  the  form 


M 


L2  it) 
0 
0 


0 

L2  (<) 
0 


0 

0 

Li  it) 


(7) 


where  we  have  simply  replaced  ki  and  k2  in.  Equation  6  with  the  dipole  decay  functions  L  i  (t)  and  L2  {t).  The 
resultant  induced  dipole  moment  for  this  definition  of  the  magnetic  polarization  dyadic  is  then 


m  (t)  =  mi  (t)  +  m2  it) 

=  Li  it)  [(z'  •  Bn  z']  +  L2  it)  [(y'  •  B'’)  y'  +  (x'  •  B^)  x'] 

Therefore,  our  approximate  forward  model  represents  the  TEM  response  of  two  orthogonal  dipoles.  The  first 
dipole  is  parallel  to  the  symmetry  axis  of  the  target,  and  the  second  dipole  is  perpendicular  to  the  symmetry  axis. 
These  dipoles  decay  independently  according  to  the  decay  laws  L 1  (t)  and  L2  it),  respectively. 


By  choosing  the  appropriate  parameters,  this  ’two-dipole’  model  produces  TEM  responses  that  are  consistent  with 
those  observed  field  measurements  of  UXO.  It  has  been  noted  that  the  shape  anomaly  of  the  measured  response 
for  UXO  changes  with  time  (Grimm,  Blohm,  and  Lavely  1997).  The  physical  phenomenon  that  gave  rise  to  the 
temporal  changes  in  shape  anomaly  was  explained  in  terms  of  the  nature  of  the  induced  eddy  currents.  Eddy 
currents  that  circulate  end-to-end  in  the  UXO  dominate  at  early  time  but  decay  away  quicMy,  while  eddy  currents 
that  circulate  about  the  long  axis  extend  later  into  time.  This  observed  field  behavior  can  be  dupficated  by  letting 
the  two  orthogonal  dipoles  mi  (t)  and  m2  (t)  decay  independently  of  each  other.  The  dipole  m  1  (t)  is  parallel  to 
the  long  axis  and  it  simulates  the  magnetic  fields  that  arise  from  currents  circulating  about  the  axis.  The  dipole 
m2  it)  is  perpendicular  to  the  long  axis  and  it  simulates  the  magnetic  fields  that  arise  from  currents  circulating 
end-to-end.  By  assigning  a  different  decay  characteristic  (governed  by  its  decay  parameters)  to  each  dipole,  the 
relative  contribution  by  each  dipole  to  the  secondary  field  can  vary  with  time. 


Time  Decay  Functions  Li(^)  and  L2{t) 

The  time  decay  for  a  sphere  is  determined  by  the  sum  of  exponentials.  This  result  generalizes  to  the  case  of 
a  conductive  body  of  arbitrary  size  and  shape  in  an  insulating  medium  illuminated  by  a  step-off  primary  field 
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(Kaufman  1994).  Thus,  the  form  for  L  (t)  should,  at  least,  be  able  to  duplicate  the  time  decay  features  observed 
for  the  sphere.  Plots  of  the  B-field  and  9B /dt  response  for  both  a  magnetically  permeable  (e.g.,  steel)  and 
nonpermeable  (e.g.,  aluminum)  sphere  are  shown  in  Figure  2. 

An  appropriate  form  of  the  decay  law  for  the  B-field  is 

L  {t)  =  k  {t  +  a)~^  e~ .  (9) 

The  parameter  k  controls  the  magnitude  of  the  modelled  response.  The  three  parameters  a,  j5,  and  7  control  the 
duration  and  characteristics  of  the  three  different  stages  of  the  time  decay  curve.  The  duration  of  the  relatively 
flat  early  time  stage  is  proportional  to  the  parameter  a.  The  linear  decrease  of  response  observed  during  the 
intermediate  time  stage  is  determined  by  t  The  exponential  decay  characterizing  the  late  time  stage  is  controlled 
by  the  parameter  7.  Figure  2a  demonstrates  the  ability  to  reproduce  the  secondary  B-field.  This  form  of  the  decay 
law,  with  the  a  parameter  absent,  was  suggested  to  us  in  a  personal  communication  fi’om  J.D.  McNeill. 


The  time  derivative  9B/5t,  measured  directly  with  most  TEM  receivers,  can  also  be  modelled  with  Equation  9. 
Figure  2a  includes  plots  of  the  9B/9t  curves  for  a  steel  and  an  aluminum  sphere.  The  early  time  behaviour 
for  the  nonpermeable  sphere  follows  a  t  decay,  and  these  curves  are  different  from  those  of  B  in  Figure  2b. 
Nevertheless,  the  curves  are  still  represented  by  early  time  tum-overs,  and  linear  and  exponential  decays  that  can  be 
accommodated  by  Equation  9.  The  suitability  is  demonstrated  by  the  fit  between  the  laboratory  measured  response 
and  a  predicted  response  obtained  by  evaluating  Equation  9. 


In  the  following  section,  we  genetically  denote  the  TEM  response  as  £  (r,  t)  where  ^  can  be  the  magnetic  field  or 
its  time  derivative.  The  time  dependent  decay  of  ^  is  given  by  Equation  9. 


The  Approximate  Forward  Model 


With  the  above  background,  we  can  write  an  approximate  expression  for  the  secondary  field  response  of  an  axi- 
symmetric  target.  First,  let  us  switch  from  the  body-fixed  (primed)  coordinate  system  to  a  space-fixed  coordinate 
system,  which  is  more  amenable  to  the  definitions  of  target  and  sensor  location  of  a  typical  field  survey  (Figure  3). 
A  vector  v'  in  the  body-fixed  coordinate  system  is  related  to  a  vector  v  in  the  space-fixed  coordinate  system  via 
the  Euler  rotation  tensor  A  9,  <f))  by  (Arfken  1985): 

v'  =  Av  (10) 


Because  of  the  axial  symmetry  of  the  problem  tp  —  0,  the  Euler  rotation  tensor  can  be  written: 


A  = 


cos  6  cos  (p 
—  sin^ 
sin  9  cos  (p 


cos  0  sin  ^  —  sin  9 

coscp  0 

sin  9  sin  p  cos  9 


(11) 


where  9  is  the  angle  between  the  symmetry  axis  of  the  target  (z'  in  Figure  lb)  and  the  vertical  axis  in  the  space- 
fixed  coordinate  system  (z  in  Figure  3),  and  (f>  is  the  angle  between  the  projection  of  z'  onto  the  horizontal  plane 
and  X. 


The  approximate  forward  modelling  is  written  by  substituting  the  definition  of  the  induced  dipole  of  Equation  8 
into  the  expression  for  a  dipole  field  (Equation  1),  and  carrying  out  the  dyadic  product.  Let  us  consider  a  target 
whose  center  is  located  at  R  in  the  space-fixed  coordinate  system.  The  secondary  response  ^  (r,  t),  measured  at 
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a  receiver/transmitter  location  r  and  at  a  time  t  after  the  termination  of  the  primary  field,  is  then  the  sum  of  the 
responses  of  the  two  orthogonal  dipoles: 


(12) 


where 


and 


Ci  (r,i) 


Mo 

Att 


^3  [mi  (<)  •  (r  -  R)] 


(r-R) 

\r--Rf 


mi{t)  ^ 

|r-R|V 


mi  (t)  =Li  {t)  (z'  •  RP)  z' 

ma  (i)  =1/2  (t)  [(x'  ■  B*")  x'  +  (y'  ■  B**)  y'] 


(13) 


(14) 

(15) 


are  the  dipole  parallel  and  perpendicular  to  the  axis  of  symmetry.  The  unit  vectors  are  given  by  Equation  11. 


In  summary,  the  approximate  response  of  a  buried  metallic  object  given  by  Equation  12  can  be  generated  from  13 
parameters  that  describe  the  object.  These  model  parameters  are  elements  of  the  model  vector 

m  =  [A",  Y,  Z,  (/),  6,  ki,  qi,  7i,A;2,  aa,  ^2,  72]  (16) 

where  X  and  F  denote  the  surface  projection  of  the  centroid  of  the  body,  and  Z  is  the  depth  of  the  object  below  the 
surface.  The  orientation  of  the  target  is  described  by  the  two  angles  9  and  (p.  The  remaining  parameters  describe 
the  decay  characteristics  of  the  two  dipoles:  A;  1 ,  ai ,  ySi ,  and  71  describe  the  dipole  parallel  to  the  axis  of  symmetry 
(mi),  and  k2,  02,  P2,  and  72  describe  the  dipole  perpendicular  to  the  axis  of  symmetry  (m2).  Thus,  the  inversion 
for  the  model  m  will  immediately  give  estimates  of  target  location  and  orientation.  Information  on  the  shape,  size, 
and  material  parameters  of  the  target  may  later  be  inferred  firom  the  remaining  parameters. 
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3  Nonlinear  Parameter  Estimation  Procedure 


In  this  paper,  it  is  first  assumed  that  the  response  measured  in  a  survey  is  the  result  of  a  single  body,  and  second, 
that  the  response  of  this  single  body  can  be  accurately  modelled  with  Equation  12.  With  these  hypotheses,  an 
inversion  procedure  is  developed  that  utilizes  the  approximate  forward  model. 


The  forward  model  can  be  expressed  as: 

dj  =  Fj[m]  and  j  — 1,2,3,. ..N  (17) 

This  equation  expresses  the  mapping  of  the  model  vector  m  to  a  datum  djhya  functional  Fj .  The  forward  mapping 
Fj  is  defined  by  Equation  12  and  is  a  nonlinear  functional  of  the  13  model  parameters  given  in  Equation  16.  In  the 
inverse  problem,  these  parameters  are  retrieved  from  a  vector  of  observed  data  d  by  minimizing  a  least-squares 
objective  function.  Before  proceeding  to  the  details  of  the  inversion,  there  are  two  important  practical  aspects  to  be 
introduced:  (1)  Selected  parameters  must  remain  positive  and  (2)  Parameters  should  be  scaled  to  enhance  stability 
in  the  iterative  process. 


In  the  approximate  forward  model,  the  time  decay  parameters  ki,ai,Pi,  and  7^  (i  =  1 , 2)  are  defined  as  positive.  In 
the  inverse  problem  the  positivity  of  these  parameters  can  be  maintained  by  solving  the  associated  square-variable 
unconstrained  problem  (Gill,  Murray,  and  Wright  1981).  Following  this  formulation,  each  time  decay  parameter 
mi  is  replaced  by  the  squared  variable  Wi  such  that  mi  =  wf.  A  second  transformation  is  a  linear  scahng  that 
ensures  that  each  component  of  m  is  of  order  unity.  This  is  done  by  dividing  each  parameter  m  ,  by  its  typical  value. 
In  the  following  presentation  of  the  inversion  algorithm,  the  letter  m  continues  to  denote  the  parameter  vector  to 
be  inverted  for,  but  now  includes  the  positivity  and  scaling  transformations. 


Defining  the  Objective  Function 


If  there  are  L  time  chaimels  and  K  locations  where  TEM  data  are  collected,  then  there  will  be  =  KL  data 
points  contained  in  the  data  vector  d  Because  data  wiU  be  collected  on  several  lines,  with  a  number  of  stations 
per  line,  there  will  generally  be  far  more  data  than  model  parameters  {N  »  13).  Therefore,  the  inversion  for  m 
involves  solving  an  overdetermined  system  of  nonlinear  equations,  with  the  goal  of  finding  the  model  that  produces 
the  data  that  best  fit  the  observed  data.  This  is  a  nonlinear  least  squares  problem  and  is  solved  by  minimizing; 

=  (18) 


where  F[m]  are  the  forward  modelled  data,  d  are  the  observed  data,  and  $  is  the  least-squares  objective  function 
that  measures  how  closely  our  predieted  data  match  the  observed  data.  W  ^  is  the  data  weighting  matrix.  If 
the  data  are  contaminated  with  unbiased  Gaussian  random  noise,  then  W  d  is  ideally  a  diagonal  matrix  whose 
elements  are  the  reciprocals  of  the  standard  deviation  of  each  datum.  The  noise  arises  from  many  sources,  including 
sensor  location  errors,  instrument  noise,  and  inaccuracy  of  the  forward  modelling.  It  is  likely  that  the  Gaussian 
independent  assumption  is  not  valid,  but  it  is  essential  to  estimate  a  quantity  that  reflects  the  uncertainty  in  each 
datum.  We  assume  that  the  errors  can  be  characterized  by  a  percentage  of  the  datum  value  plus  a  threshold,  that  is 


i^d)ii  = 


1 

pdi  -t-  e 


(19) 
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where  p  is  typically  a  percentage  and  e  is  a  constant  that  characterizes  ambient  noise.  The  positive  e  ensures 
that  small  data  points  would  have  reasonable  errors  assigned  to  them,  and  thus  prevents  them  from  having  undue 
influence  on  the  solution. 


Minimizing  the  Objective  Function 


For  ease  of  notation,  the  least-squares  problem  is  rewritten  as 


minimize  $  (m) 


i=l 


(20) 


where  r  is  the  residual  function 

r(m)=Wd  (F[m]-d°'’*) 

andri  (m)  is  the  i*'*  component  of  r  (m).  We  adopt  a  modified  Newton’s  method  to  minimize  the  objective  function. 
The  approach  taken  here  is  to  first  make  an  initial  guess  of  the  model  parameters  m  Techniques  for  making  this 
initial  guess  are  outhned  in  Pasion  (1999).  The  starting  model  ruo  is  iteratively  improved  to  find  the  minimum 
of  Equation  20.  At  each  iteration  a  Newton’s  search  direction  is  chosen  that  minimizes  the  local  quadratic  model 
about  the  current  iterate  m  * .  The  Newton  step  <5m  for  the  nonlinear  least-squares  problem  is  then  given  by: 

H(mfc)5m= -J(m*)^r(mfc)  (21) 

where  the  Jacobian  matrix  J  and  the  Hessian  matrix  H  are  defined  as 

=  and  H  =  J  (m)^  J  (m) -f  S  (m) 

arrij 

where  S  (m)  is 


N 

S(m)  =  ^ri(m)V^ri(m) 

i=:l 

Anewmodel  m*+i  =  m/t -I- A(im  is  then  defined,  where  the  positive  scalar  A  is  chosen  such  that  #  (m*  -f  A(5m)  < 
$  (m*).  The  sequence  of  iterations  is  terminated  once  the  relative  gradient  measure  is  less  than  a  tolerance  level , 
or  once  there  is  insignificant  change  in  the  models  obtained  at  successive  iterations  (Dennis  and  Schnabel  1983). 

Error  Bounds  of  the  Parameter  Estimates 

Once  the  model  parameters  m  *  which  minimize  the  objective  function  #  (m)  have  been  obtained,  we  can  examine 
the  reliability  and  precision  of  the  estimated  parameters  via  the  model  covariance  matrix.  Let  m  be  the  best 
estimate  of  the  model  in  the  absence  of  noise,  and  <5m  ,  =  m+  —  m«.  The  model  covariance  matrix  is  defined 
as  the  expectation  value  of  (im  *  (im,^(Bard  1974): 

Vm  =  E  (Jm*  (im  J)  E  (h*-i  J  J 

The  Hessian  and  Jacobian  in  the  above  expression  are  evaluated  at  m  =  m*,  and  are  therefore  constants.  As  a 
result  they  can  be  taken  outside  of  the  expectation  value  expression: 

V„  =  H*-'J/ Vd  (22) 
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where  Vd  is  the  covariance  matrix  of  the  data.  In  the  case  when  the  observations  have  uncorrelated  errors,  the 
data  covariance  matrix  reduces  to  a  diagonal  matrix,  and  an  estimate  of  the  standard  deviation  of  the  i  model 
parameter  is  then: 

OT  =  =-  fn.-'  J  (23) 

V  /  a 

Model  variance  estimates  applied  to  nonlinear  problems  are  not  as  reliable  as  when  implemented  in  hnear  least- 
squares  problems,  and  they  should  only  be  used  as  a  very  rough  estimate  (Bard  1974;  Dennis  and  Schnabel  1983). 
Nevertheless,  Equation  23  at  least  provides  a  minimum  estimate  to  the  uncertainties  of  the  parameters. 
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4  Relating  Model  Parameters  to  Material  and 
Geometric  Properties 


The  above  inversion  generates  the  parameters  that  characterize  a  target’s  TEM  anomaly.  The  next  step  is  to  interpret 
these  parameters.  Recall  that  UXO  are  typically  rod-like  rather  than  plate-like  and  are  magnetically  permeable. 
In  order  to  extract  these  potentially  UXO  identifying  features  from  the  recovered  model  m  *,  we  use  the  inversion 
procedure  to  fit  a  series  of  decay  curves  from  a  range  of  axi-symmetric  targets  of  different  shape,  geometry,  and 
material  properties.  We  then  generate  empirical  relationships  between  the  parameters  and  target  characteristics. 
The  data  curves  used  for  this  analysis  were  either  TEM  measurements  made  in  the  Geonics  Ltd.  laboratory,  or  they 
were  synthetically  generated  decay  curves  for  a  sphere  using  Equation  12. 


Lab  Setup  and  Measurements 

A  series  of  TEM  measurements  of  metallic  targets  was  made  by  Geonics  Ltd.  A  40  m  x  40  m  square  transmitter 
loop  was  used  to  provide  a  relatively  uniform  field  at  the  center  of  the  loop.  Aim  diameter  receiver  coil  was 
placed  coaxial  and  coplanar  to  the  transmitter  loop,  and  each  target  was  located  at  the  center  of  the  receiver  loop. 
The  Geonics  PROTEM  47  time  domain  equipment  was  used  for  producing  the  transmitting  field  and  for  recording 
the  time  domain  measurement  as  a  result  of  a  step-off  current.  Measurements  of  the  time  decay  response  of  these 
targets  were  recorded  as  plots  of  log(9B /dt)  vs  log(t).  Since  values  were  not  recorded  by  a  data  logger,  the  plots 
were  subsequently  digitized.  Plots  of  the  steel  target  responses  were  digitized  by  hand  by  J.D.  McNeill  Geonics 
Ltd.,  and  the  aluminum  target  responses  were  digitized  at  UBC  after  scanning  the  plots  into  a  computCT.  Analyses 
were  performed  on  both  the  impulse  (5B  /  dt)  response  measured  as  an  induced  voltage  in  the  receiver  and  also  on 
the  B-field  response.  The  B-field  response  was  obtained  by  integrating  the  induced  voltage. 


Two  sets  of  targets  were  measured.  The  first  set  of  measurements  involved  recording  the  TEM  response  for  a  series 
of  steel  and  aluminum  rectangular  prisms  of  different  dimensions.  Each  prism  had  at  least  one  dimension  of  20.3 
cm,  and  the  targets  ranged  from  a  thin  rod  (20.3  x  0.6  x  0.6  cm)  to  a  cube  (20.3  x  20.3  x  20.3  cm)  to  a  thin 
plate  (20.3  x  20.3  x  0.6  cm).  A  second  set  of  measurements  was  made  on  24  sample  UXO.  These  targets  included 
various  ordnance  items  used  by  NATO  since  World  War  n.  The  ordnance  ranged  in  length  from  18  to  85  cm  and 
in  diameter  from  6.1  to  15.9  cm.  A  diagram  of  all  the  ordnance,  along  with  a  table  listing  the  dimensions  of  each 
ordnance,  is  included  in  Pasion  (1999). 


The  axi-symmetric  targets  were  placed  in  two  orientations  at  the  center  of  the  receiver  loop.  Each  target  was 
measured  with  the  axis  of  symmetry  perpendicular  and  parallel  to  the  primary  field.  Since  the  strength  of  each 
induced  dipole  is  proportional  to  the  projection  of  the  primary  field  onto  the  dipole  direction,  the  two  measurement 
orientations  isolate  the  decay  behaviour  of  each  of  the  two  dipoles.  For  example,  consider  a  plate.  When  the 
primary  field  is  perpendicular  to  the  plane  of  the  plate,  the  projection  of  the  primary  field  onto  dipole  2  is  zero, 
thus  the  approximate  forward  model  assumes  the  response  can  be  modelled  as  a  single  dipole  perpendicular  to  the 
plate.  The  decay  parameters  of  dipole  1  (fc  i ,  ai ,  /3i ,  and  71 )  can  then  be  estimated  by  fitting  this  curve  to  the  decay 
law  (Equation  9).  When  the  primary  field  is  parallel  to  the  plane  of  the  plate,  the  response  is  the  result  of  the  fact 
that  dipole  2  and  parameters  k2,a2,  and  72  can  be  recovered. 
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Relationships  between  the  target  characteristics  and  the  model  parameters  were  established  in  the  following  man¬ 
ner.  A  scaled-down  version  of  the  nonlinear  least-squares  techniques  outlined  in  the  previous  section  was  used 
to  obtain  the  decay  parameters  k,  a,  p,  and  7  for  each  of  the  target’s  two  dipoles.  Secondly,  we  observed  how 
recovered  values  of  model  parameters,  or  combinations  of  parameters,  changed  with  the  dimensions  and  mag¬ 
netic  properties  of  the  measured  prism.  The  patterns  in  the  behaviour  of  the  parameters  then  led  to  the  shape  and 
permeability  discrimination  diagnostics  that  are  proposed  in  the  following  sections. 


A  Relationship  Between  /3  and  Magnetic  Permeability 

UXO  are  generally  made  of  steel,  which  is  a  ferrous  material.  Therefore,  the  magnetic  permeability  is  likely  to 
be  an  identifying  characteristic  of  UXO.  To  generate  a  link  between  magnetic  permeability  and  model  parameters, 
forward  modelled  responses  were  calculated  for  a  series  of  spheres  varying  in  size  and  permeability.  Both  B-field 
and  data  were  then  inverted  to  generate  decay  parameters,  and  in  particular,  to  produce  estimates  of  the 

parameter  /?.  The  plots  of  as  a  function  of  sphere  radii  and  magnetic  permeability  are  provided  in  Figure  4. 


Figure  4  suggests  that  the  value  of  /3  obtained  for  a  sphere  may  be  diagnostic  in  determining  whether  the  sphere 
is  permeable  or  nonpermeable.  Figure  4a  exhibits  the  relationship  for  the  5B /dt  responses.  For  a  sphere  with 
fir  =  150  (typical  for  steel),  we  see,  for  spheres  with  radii  between  5  to  15  cm,  that  /?  falls  between  1.11  and  1.35; 
while  for  a  sphere  with  fir  =  1,  P  has  a  value  of  approximately  0.5,  which  corresponds  to  the  early  time  t 
behavior  that  Kaufman  (1994)  predicted  for  a  nonpermeable  sphere.  Therefore,  when  applying  our  inversion  to 
the  time  derivative  of  the  field,  a  value  of  greater  than  about  0.8  indicates  that  the  target  is  most  likely 

permeable.  This  analysis  is  repeated  on  the  forward  modelled  B-field  responses,  and  the  results  are  plotted  in 
Figure  4b.  A  threshold  value  of  P^  =  0.3  could  be  used  such  that  targets  with  a  B-field  response  characterized  by 
a  /3  >  0.3  indicates  a  permeable  target. 


The  use  of  /9  as  a  diagnostic  to  determine  permeability  can  be  extended  to  nonspherical  targets  by  looking  at  the 
recovered  p  values  for  the  aluminum  and  steel  prisms.  The  inversion  produces  two  values  of  /?,  one  for  each  of  the 
excited  dipoles,  to  describe  a  buried  target.  We  suggest  taking  the  average  of  the  two  recovered  P  values,  which  we 
label  as  p.  When  analyzing  the  9B/9t  responses  of  the  axi-symmetric  aluminum  targets,  =  0.52  with  a 

standard  deviation  of  0.07.  For  the  steel  targets,  ^  standard  deviation  of  0.08.  These  averages 

fall  on  either  side  of  the  0.8  threshold  obtained  by  fitting  sphere  5B /dt  responses. 

When  analyzing  the  B-field  responses  ofthe  axi-symmetric  aluminum  targets,  =  0.1 7  with  a  standard  deviation 
of  0.03.  For  the  steel  targets  P^  =  0.5  with  a  standard  deviation  of  0.2.  These  averages  fall  on  either  side  of  the 
0.3  threshold  obtained  by  fitting  sphere  B-field  responses,  and  so  again,  a  consistent  criterion  can  be  used. 


Relationships  Between  Model  Parameter  Ratios  and  Target  Shape 

Empirical  relationships  were  also  established  between  the  target  shape  and  the  ratios  k  1  /fc2  and  PilP2-  For  space 
reasons,  we  present  only  the  analysis  of  the  9B /dt  response  of  the  targets,  and  refer  the  reader  to  Pasion  (1999)  for 
the  analysis  of  the  B-field  data.  There  we  show  that  the  same  target  shape  diagnostics,  developed  here  for  9B /dt 
data,  also  apply  to  B-field  data. 
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The  Ratio  ki/k2 

The  recovered  k  values  for  targets  ranging  from  a  steel  plate  to  a  steel  rod  are  shown  in  Figure  5a,  and  the  calculated 
jfc-ratios  are  shown  in  Figure  5b.  For  a  steel  plate,  the  fc-ratio  kijk^  <  1.  For  a  steel  bar  the  fc-ratio  ki/k-i  >  1. 
The  recovered  k  values  for  aluminum  targets  are  shown  in  Figure  5c.  The  opposite  orientation  effect  was  observed 
for  an  aluminum  rod,  i.e.,  kijk^  <  1  (Figure  5d). 


The  Ratio  A  7/52 

In  addition  to  the  relative  strength  of  the  dipoles  being  shape  dependent,  the  slope  of  the  time  decay  response  (either 
5B /dt  or  B-field)  during  the  intermediate  time  stage  is  dependent  upon  the  target  shape.  This  effect  was  seen  in 
steel  targets  only.  The  steepness  of  the  response  during  the  intermediate  time  stage  is  reflected  in  the  parameter  (3. 
The  recovered  P  values  for  targets  ranging  from  a  steel  plate  to  a  steel  rod  are  shown  in  Figure  6a,  and  /?  values  for 
aluminum  targets  are  shown  in  Figure  6c.  A  dipole  that  decays  at  a  greater  rate  will  have  a  larger  p.  The  rate  of 
decay  of  the  SB  jdt  response  is  greater  when  the  plane  of  a  steel  plate  is  perpendicular  to  the  primary  field  (dipole 
1)  than  when  the  plane  of  a  steel  plate  is  parallel  to  the  primary  field  (dipole  2).  Thus,  for  a  steel  plate,  the  /3-ratio 
/3i  7/32  >  1.  In  the  case  of  a  rod,  the  clB  /dt  response  decays  faster  (and  thus  /3  is  larger)  when  the  main  axis  of  the 
rod  is  perpendicular  to  the  primary  field  (dipole  2).  In  the  case  of  a  steel  rod,  the  /3-ratio  P1IP2  <  1  (Figure  6b). 

For  aluminum  targets  the  response  shape  looks  essentially  the  same  for  each  of  the  targets.  The  3B/ dt  response 
exhibits  a  power  law  decay  of  and  is  exponential  at  later  times.  The  decay  curves  for  alumimun  targets  are 
essentially  the  same  regardless  of  target  shape,  and  therefore  there  is  no  relationship  between  the  ^-ratio  and  the 
aspect  ratio  (Figure  6d). 


The  Discrimination  Algorithm  Using  dB/dt  Data 

The  results  from  the  previous  section  suggest  the  following  algorithm  for  using  dB/dt  data  to  help  identify  possi¬ 
ble  UXO  targets: 


1.  Perform  the  nonlinear  inversion  outlined  in  the  previous  section  to  recover  model  parameters  for  the  two- 
dipole  model. 

2.  Compute  P  =  ^  [0-^  +  P2).ll  P  >  0.80,  then  the  target  is  most  likely  permeable. 

3.  Compute  ratios  P1/P2  and  fci /fc2-  There  are  two  options: 

•  ^  >  0.8  =>  Ferrous  Target.  If  ^1/^2  >  1  and  0x102  <  1,  then  a  permeable  rod-like  target  was 
measured.  If  fci  /fc2  <  1  and  01  (02  >  1,  then  a  permeable  plate-like  target  was  measured. 

•  ^  <  0.8  =>  NonFerrous  Target  If  fci/fc2  >  1  then  nonpermeable  plate-like  target  was  measured.  If 
fci  / ^2  <  1  then  the  target  is  rod-like.  Pi  102  does  not  give  supporting,  or  extra,  information. 


The  above  algorithm  can  be  extended  to  the  analysis  of  B-field  data  simply  by  changing  the  0  threshold  to  0.3. 

The  parameters  a  and  7,  which  respectively  are  representative  of  a  target’s  early  and  late  time  characteristics,  are 
not  used  in  the  discrimination  algorithm.  The  a  parameter  is  very  sensitive  to  how  early  the  measurement  of  the 
TEM  sensor  extends  in  time,  and  we  have  found  that  it  is  not  as  robustly  estimated  as  k,  0,  and  7.  The  late  time 
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constant  for  conductors  is  related  to  a  target’s  shape,  size,  and  conductivity  (Kaufman  1994  and  Nabighian  and 
Macnae  1991).  However,  as  noted  by  Nabighian  and  Macnae  (1991),  using  the  time  constant  for  discrimination  is 
difficult  without  prior  knowledge  of  the  shape  and  either  a  length  scale  or  conductivity.  Our  tests  thus  far  have  been 
unable  to  generate  a  consistent  relationship  between  the  time  constant  recovered  from  lab  TEM  measurements  and 
target  characteristics. 
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5  Synthetic  Data  Set 


The  parameter  estimation  procedure  is  now  tested  on  a  synthetically  generated  field  data  set.  The  object  of  interest 
is  a  75  mm  antitank  mortar.  The  primary  decay  curves  for  the  axial  and  perpendicular  orientations  were  obtained 
fi:om  measurements  made  by  Geonics  using  the  setup  described  in  the  previous  section.  These  were  inverted  to 
recover  the  decay  parameters  for  each  dipole.  For  this  simulation,  the  target  is  assumed  to  be  buried  at  a  depth  of 
57  cm  (Z=l  m),  and  located  at  (2m  N,2m  E)  on  the  survey  grid.  The  mortar  is  oriented  such  that  </>  =  30  °  and 
6*  =  65  °.  Equation  12  was  used  to  generate  the  “observed”  data  for  this  example. 


The  survey  consists  of  a  2  x  2  m  grid,  containing  five  lines  of  data  separated  at  50  cm  line  spacings,  with  stations 
located  at  20  cm  intervals  along  each  line.  At  each  station,  the  vertical  component  of  the  voltage  is  generated  for 
26  logarithmically  spaced  time  channels.  The  tune  channels  range  from  0.01  to  100  ms.  In  order  to  make  this 
example  closer  to  a  real  TEM  data  set,  5  percent  random  Gaussian  noise  was  added  to  the  data  and,  since  a  real 
TEM  instrument  will  have  a  finite  measurement  sensitivity  to  the  secondary  field,  a  data  threshold  of  0.(X)1  is  set. 
There  are  1,278  total  data  points  exceeding  the  minimum  threshold. 


The  inversion  is  carried  out  with  a  data  weighting  matrix  in  Equation  19  with  p  —  0.05  and  e  =  0.001  mV.  The 
observed  Hata^  and  data  predicted  by  the  recovered  model,  are  compared  in  Figures  7  and  8.  Figure  8  shows  a  plan 
view  comparison  for  3  of  the  26  time  chaimels.  The  difference  maps  exhibit  a  random  distribution  over  the  data, 
indicating  the  reluctance  of  the  inversion  to  fit  the  noisy  portion  of  the  data.  A  comparison  of  the  true  model  m  true> 
recovered  model  mrec,  and  the  estimated  model  standard  deviations  ex’"  are  found  in  Table  1.  The  discrimination 
algorithm,  when  applied  to  the  recovered  decay  parameters,  yields  the  following.  The  value  of  /3  =  1.07(>  0.8) 
indicates  that  the  target  is  likely  permeable.  The  ratios  fc  1/^2  =  3.71(>  l)and/?i//32  =  0.71(<  1)  indicate,  for  a 
magnetically  permeable  target,  that  the  TEM  response  is  likely  from  a  rod-like  target. 
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Field  Data  Set 


We  now  apply  our  algorithm  to  a  TEM  field  data  set  acquired  at  the  ERDC  UXO  test  site  in  Vicksburg,  MS.  The 
Geonics  EM63  instrument  used  for  the  survey  is  a  multi-time  channel  time  domain  unit  consisting  of  a  1  x  1  m 
square  transmitter  coO  and  a  coaxial  horizontal  circular  receiver  loop  mounted  on  a  two-wheel  trailer.  Measured 
voltages  are  averaged  over  26  geometrically  spaced  time  gates,  spanning  the  range  0.18  to  25.14  ms. 

A  105  mm  projectile  is  placed  in  the  ground  with  its  center  at  2.0  m  East,  1 .83  m  North  and  at  a  depth  of  0.44  m 
from  the  surface.  The  projectile  was  placed  horizontal  (6  =  90  °),  with  its  tip  pointing  to  the  North  ((/>  =  0  °).  Once 
the  target  was  placed  in  the  ground,  it  was  not  covered  in  soil.  The  survey  consisted  of  a  2  x  2  m  grid  centered  on 
the  target,  containing  five  lines  running  North-South  separated  at  50  cm  line  spacing,  with  stations  located  at  5  cm 
intervals  along  each  fine.  A  measured  signal  of  less  than  1  mV  is  assumed  to  be  indistinguishable  from  the  noise. 
The  resulting  data  set  contains  1,882  total  data  points. 


The  inversion  is  carried  out  with  a  data  weighting  matrix  in  Equation  19  with  p  =  0.05  and  e  =  1  mV.  The  first 
stage  of  the  time  decay  evident  in  Figure  2  is  not  observed  in  the  time  window  recorded  by  the  EM63.  Therefore, 
we  invert  these  data  by  setting  q  to  be  a  small  constant  that  does  not  affect  the  predicted  data  within  the  EM63 
time  range.  So  only  parameters  k,  /?,  and  7  for  each  dipole  are  recovered.  The  observed  and  predicted  data  are 
compared  in  Figures  9  and  10.  Figure  10  shows  a  plan  view  comparison  for  5  of  the  26  time  channels.  At  early 
times,  the  anomaly  has  a  single  peak  located  approximately  above  the  UXO  center.  This  peak  splits  into  two 
distinct  peaks  at  late  time.  The  recovered  model  predicts  data  that  reflects  this  behavior.  Figure  9  compares  the 
decay  curve  measured  at  four  stations  on  the  survey. 


The  recovered  location  and  orientation  parameters  are  fisted  in  Table  2a.  The  recovered  easting  of  2.04  m  differs 
from  the  true  value  of  2.00  m  by  4  cm.  The  recovered  northing  of  1.77  m  differs  from  the  true  value  of  1 .83  cm 
by  6  cm,  placing  the  inducing  dipole  closer  to  the  projectile  tail.  These  errors  are  of  the  same  magnitude  as  can  be 
expected  in  spotting  the  station  location  in  the  field  survey.  In  addition,  the  buried  105  mm  projectile  has  a  copper 
rotating  band  near  the  tail  of  the  projectile.  It  has  been  suggested  that  the  presence  of  the  rotating  band  will  shift  the 
location  of  the  induced  dipole  from  the  target  center  toward  the  tail  (Miller  2000).  The  recovered  burial  depth  of 
0.47  m  is  3  cm  deeper  than  the  expected  depths  of  0.44  m.  The  orientation  parameters  9  and  0  are  well  recovered. 
The  recovered  decay  parameters  are  listed  in  Table  2b,  and  the  diagnostics  applied  to  these  parameters  are  fisted  in 
Table  2c.  The  value  of =  0.91(>  0.8)  indicates  that  the  target  is  Ukely  to  be  magnetically  permeable.  The  ratios 
k\/k2  =  2.63(>  1)  and/)i//?2  =  0.69(<  1)  indicate,  for  a  magnetically  permeable  target,  that  the  TEM  response 
is  likely  to  be  from  a  rod-like  target. 
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7  Discussion  and  Conclusion 


Efficient  remediation  of  areas  containing  UXO  first  requires  that  purely  conductive  metal  targets  be  distinguished 
from  steel  targets  that  are  conductive  and  permeable.  The  second  stage  then  focuses  upon  determining  if  the 
steel  target  is  rod-like  (and  likely  to  be  a  UXO)  or  plate-like  (and  not  be  of  interest).  To  attack  this  problem,  we 
propose  a  modified  parametric  model  from  which  TEM  responses  can  be  estimated.  The  TEM  response  of  a  buried 
axisymmetric  metallic  object  is  modelled  as  the  sum  of  two  dipoles  located  at  the  midpoint  of  the  body.  Nonlinear 
inversion  methods  are  used  to  extract  the  parameters  from  the  field  data  and  these  parameters  are  subsequently 
used  in  a  discrimination  procedure  that  has  two  parts.  First,  the  decision  about  whether  the  object  is  ferrous  might 
be  made  by  examining  the  size  of  the  recovered  P  values.  Second,  if  the  object  is  considered  to  be  ferrous,  then 
the  ratios  of  A:i/A:2  and  P1IP2  are  diagnostic  indicators  of  whether  the  geometry  is  plate-Uke  or  rod-like.  These 
diagnostics  were  developed  for  both  B -field  and  5B  /  dt  TEM  data. 


This  algorithm  was  applied  to  a  synthetic  data  set  as  well  as  to  a  field  data  set  collected  by  the  Geonics  EM63  time 
domain  electromagnetic  sensor  over  a  105  mm  projectile.  In  both  cases  the  diagnostics,  applied  to  the  recovered 
model  parameters,  correctly  predicted  that  the  TEM  anomaly  was  produced  by  a  magnetically  permeable  and  rod- 
like.  metallic  target.  Although  further  testing  will  be  required  to  fuUy  evaluate  our  proposed  technique,  the  results 
presented  here  are  promising  and  may  have  a  positive  impact  on  the  interpretation  of  UXO  detection  data. 
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Figure  1:  The  body-fixed  (primed)  coordinate  system  for  a  sphere  and  a  spheroid. 


(a)  The  time  decay  behaviour  of  the  time  derivative  of  the 
magnetic  field  cJB  /  dt. 


(b)  The  time  decay  behaviour  of  the  magnetic  flux  density 
B. 


Figure  2:  Fitting  the  time  decay  of  a  sphere  using  L{t).  The  B-field  response  is  normalized  by  the  strength  of  the  primary  field. 
The  solid  lines  are  responses  evaluated  from  eq.  (9).  The  agreement  supports  the  validity  of  this  parametric  representation  of 
the  time  domain  responses. 
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Figure  3;  The  field  (unprimed)  coordinate  system  for  a  buried  target.  The  unit  vectors  x,  y,  and  z  define  the  field  co-ordinate 
system,  and  x',  y',  and  z'  define  the  body-fixed  coordinate  system. 


(a)  Recovering  f)  by  fitting  the  d^/dt  data. 


(b)  Recovering  /?  by  fitting  the  B-field. 


Figure  4:  The  behaviour  of  parameter  j3  for  various  size  spheres  with  var5dng  permeabibty  fi. 


(a)  Recovered  k  parameter  from  fitting  the  mea-  (b)  Relationship  between  the  ki/k2  ratio  derived  from  /dt  data  and 

suied  9B  /  dt  response  of  steel  axi-symmetric  tar-  the  shape  of  a  steel  target 

gets. 


(c)  Recovered  k  parameter  fiom  fitting  the  mea-  (d)  Relationship  between  the  fcj  / k-z  ratio  and  the  shape  of  an  aluminum 

sured  d'B/dt  response  of  aluminum  axi-symmetric  target. 

targets 


Figure  5:  Relating  the  aspect  ratio  of  a  steel  target  with  the  ratio  ki /k2 
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8  X  1/4  X  1/4  I - 1  I  I - 1 - 1 - 1 - r8x1/4x1/4 


(a)  Recovered  /3  parameter  firom  fitting  the  mea-  (•>)  Relationship  between  the  ft  ratio  derived  firom  d^/dt  data  and 

sured  jdt  response  of  steel  axi-symmetric  tar-  shape  of  a  steel  target. 

gets. 


(c)  Recovered  ^  parameter  from  fitting  the  mea-  (d)  Relationship  between  the  ft  //S2  ratio  and  the  shape  of  an  aluminum 

sured  9B/5t  response  of  aluminum  axi-symmetric  target. 

targets. 


Figure  6:  Relating  the  aspect  ratio  of  a  steel  target  with  the  ratio  /3i  Ifi2 
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exIMxIM  I - 1 - 1  8x174x1/4 


Table  1:  Recovered  parameters  from  the  inversion  of  the  synthetic  data  set.  mo  is  the  starting  model.  The  true  model  mtme  and 
the  recovered  model  mrec  are  close.  The  difference  between  the  true  and  recovered  model  falls  within  the  estimated  standard 
deviation. 

a.  Location  and  Orientation  b.  Decay  Parameters 


nrii 

mo 

^true 

rrirec 

a"' 

Northing  (m) 

1.90 

2.00 

2.00060 

0.00254 

Easting  (m) 

2.15 

2.00 

2.00105 

0.00169 

Depth  from 
loop  (m) 

1.00 

1.00051 

0.00764 

4>  (degrees) 

30 

30.07 

0.18 

6  (degrees) 

45 

65 

65.04 

0.11 

rrii 

m<, 

^true 

^rec 

h 

7.07 

12.02 

12.064 

0.18 

ai 

0.01 

0.0076 

0.00759 

0.0030 

/3i 

1.00 

0.89 

0.890 

0.0070 

7i 

17.635 

ifJHM 

0.071 

Oto 

0.01 

0.0077 

0.0076 

0.0027 

1^2 

1.00 

1.25 

1.252 

0.014 

72 

3.16 

11.54 

11.68 

0.33 

Voltage 

(mV) 


Figure  7;  The  observed  and  predicted  decay  curves  for  four  stations  in  the  synthetic  data  set  inversion.  The  predicted  vertical 
component  of  the  response,  represented  by  the  solid  lines,  are  a  good  match  to  the  artificially  generated  noisy  data  set. 
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Figure  9:  The  observed  and  predicted  decay  curves  for  four  stations  in  the  105  mm  projectile  UXO  field  data  set  inversion. 
The  predicted  decay  of  the  vertical  component  of  the  measured  voltages  are  represented  by  the  solid  lines,  and  the  symbols 
represent  the  Geonics  EM63  field  measurements. 


26 


Chapter  7  Discussion  and  Conclusion 


Table  2:  Recovered  parameters  for  the  field  data  inversion,  a.  demonstrates  that  the  inversion  was  successful  in  obtaining  the 
approximate  location  and  orientation  of  the  target,  b.  lists  the  recovered  decay  parameters  of  the  two  dipoles,  c.  lists  the  results 
of  applying  the  identification  diagnostics  to  the  recovered  decay  parameters.  Application  of  the  diagnostics  indicates  that  the 
buried  target  is  permeable  and  rod-like  and  therefore  a  candidate  for  UXO. 

a.  Location  and  Orientation  b.  Decay  Parameters 


mo 

^rec 

Expected 

Parameters 

a"’ 

Northing  (m) 

1.9 

1.77 

1.83 

0.008 

Easting  (m) 

2.15 

2.04 

2.00 

0.005 

Burial 
Depth  (m) 

0.6 

0.47 

0.44 

0.01 

4>  (degrees) 

45 

10.1 

~  0 

0.9 

d  (degrees) 

45 

84.7 

~90 

0.13 

m 

mo 

cr'" 

ki 

7.07 

76.8 

2.5 

Pi 

1.00 

0.74 

0.04 

7i 

3.16 

31.8 

8.3 

k2 

7.07 

29.2 

3.7 

P2 

1.00 

1.08 

0.07 

72 

3.16 

6.1 

1.3 

c.  Diagnostics 


Diagnostic 

Result 

Conclusion 

P 

0.91 

permeable 

ki/k2 

2.63 

rod-like 

P1/P2 

0.69 

rod-like 
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Part  II 


Verification  of  the  Two-Dipole  Modelling  and 
Inversion  Algorithm 


1  EM63  Data  Collection  at  the  U.S.  Army 
Engineer  Research  and  Development  Center 
(ERDC)  Unexploded  Ordnance  (UXO)  Test  Site 


Time  domain  electromagnetic  (TEM)  measurements  were  carried  out  between  April  5  and  April  21,  2000,  at  the 
U.S.  Army  Engineering  Research  Development  Centre  (ERDC),  Vicksburg,  MS,  to  evaluate  the  apphcability  and 
performance  of  the  forward  modelling  and  inversion  algorithm  outlined  in  Part  I  of  this  report.  The  Geonics  EM63 
TEM  sensor  was  used  for  all  the  data  collection  in  this  investigation.  The  EM63  is  a  multi-time  channel  time 
domain  unit  consisting  of  a  1  x  1  m  square  transmitter  coil  and  a  co-axial,  horizontal  50  cm  diameter  receiver 
coil  mounted  on  a  two-wheel  trailer.  More  information  about  the  EM63  sensor  can  be  found  in  “EM63  Full  Time 
Domain  Electromagnetic  UXO  Detector:  Operating  Instructions”  (2000). 


Plan  Measurements 

A  series  of  surveys  were  conducted  over  single  targets  seeded  at  the  ERDC  UXO  test  site  (Figure  1 1).  A  4  x  4  m 
square  (centered  at  18  m  A,  24  m  £)  was  chosen  for  the  surveys.  Prior  to  seeding  the  individual  targets,  EM63 
and  EM61-HH  surveys  were  carried  out  to  ensure  that  the  area  was  “quiet”,  i.e.  to  ensure  the  area  did  not  contain 
metallic  scrap.  The  borders  of  the  4  x  4  m  square  were  marked  with  string  to  indicate  the  extent  of  the  survey  area. 
Lines  for  each  survey  were  run  in  a  A  to  S  direction  with  a  line  spacing  of  50  cm.  The  location  of  the  sensor  was 
measured  more  accurately  by  marking  surv^  lines  at  1  m  spacing  with  string  and  by  dropping  a  plumb  line  from 
the  center  of  the  receiver/transmitter  loop  pair.  The  EM63  was  set  to  record  a  time  decay  curve  at  10  cm  intervals 
triggered  by  the  odometer  in  the  EM63  trailer  wheels. 


Targets  measured  in  this  survey  setup  included  several  UXO  (37  mm  projectile,  60  mm  mortar,  81  mm  mortar, 
stokes  mortar,  105  mm  mortar  projectile,  155  mm  projectile)  and  a  variety  of  scrap  items  excavated  during  a  UXO 
remediation  project  at  Camp  Croft,  Maryland.  These  targets  were  placed  at  approximately  the  center  of  the  grid 
(2  m  A,  2  m  E),  and  at  depths  up  to  75  cm.  In  all  cases,  the  strike  of  the  target  was  parallel  to  the  line  direction. 
In  order  to  save  time,  the  soil,  removed  when  digging  a  hole  for  the  target,  was  not  replaced  over  the  target.  A 
wooden  plank  was  placed  over  the  hole  in  which  the  object  was  laid.  Targets  were  generally  measured  in  three 
orientations:  horizontal,  vertical,  and  an  intermediate  angle. 


Decay  Measurements 

Controlled  experiments  were  carried  out  to  carefully  examine  how  the  secondary  field  of  a  target  decays  as  a 
function  of  depth  and  orientation.  For  these  experiments,  measurements  were  required  with  the  EM63  transmit¬ 
ter/receiver  cod  directly  above  the  center  of  various  targets  positioned  at  several  depths  and  orientations. 


To  accurately  and  quickly  position  each  target,  Mr.  Jose  Llopis,  ERDC,  designed  and  budt  a  target  holder  (Fig- 
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ure  12a).  The  jig  was  made  of  wood  and  glue,  contained  no  metallic  materials,  and  could  orient  each  target  at  15 
deg  increments  from  vertical  to  horizontal.  The  size  of  the  jig  made  it  difficult  to  bury.  Therefore  the  jig  was  only 
partially  buried  and  2x6  in.  planks  were  used  to  adjust  the  height  between  the  EM63  and  the  target  (Figure  12b). 
The  height  of  the  planks  was  varied  by  changing  the  number  of  cinder  blocks  used  to  elevate  the  planks.  The 
majority  of  measurements  were  taken  at  two  sensor  heights:  Z  ~  50  cm  and  Z  ~  1(K)  cm  from  the  center  of  the 
receiver  to  the  center  of  the  target,  where  Z  is  the  vertical  distance  from  the  center  of  the  target  to  the  center  of  the 
receiver  loop.  Additional  measurements  with  Z  ~  75  cm  were  made  on  a  subset  of  the  targets.  TEM  soundings 
were  recorded  for  several  UXO  items  ranging  from  a  37  mm  to  a  105  mm  projectile,  as  well  as  several  scrap  items. 
Samples  soundings  are  plotted  in  Figure  13. 


Several  of  the  smaller,  non-UXO  items  did  not  produce  a  significant  response  when  placed  in  the  jig.  Therefore, 
they  were  measured  using  the  setup  of  Figure  14a.  Each  target  was  placed  in  two  orientations  and  at  ground 
level.  To  facilitate  the  collection  of  cleaner  data,  a  pair  of  2  x  6  in.  planks  was  placed  on  the  ground  to  provide 
a  level  surface  for  the  EM63  to  be  pulled  along.  The  measured  voltage  curves  for  scrap  targets  1  to  8  are  plotted 
in  Figures  14b  and  c.  The  line  profiles  for  the  first  time  channel  (t  =  180  [is)  and  the  tenth  time  channel  (t  = 
0.72  msec)  are  plotted  in  Figure  15. 

EM63  TEM  Survey  of  the  ERDC  UXO  Test  Site 

The  ERDC  UXO  Test  Site  eonsists  of  a  30  m  x  95  m  grassy  field  containing  seeded  UXO  and  scrap.  On  the 
eastern  end  of  the  site  there  are  approximately  25  targets  buried  along  four  north-south  oriented  rows.  An  EM63 
surv^  was  carried  out  between  75  m  and  93.5  m  E.  The  survey  was  completed  with  1  m  line  spacing  and  10  cm 
station  spacing.  To  maintain  straight  lines,  string  was  laid  at  2  m  intervals.  Two  data  sets  were  obtained.  The 
first  data  set  was  from  75  m  to  86.5  m  E,  and  the  second  data  set  was  firom  82  m  E  to  93.5  m  E.  Figures  16a 
and  b  contain  contour  plots  for  the  first  (0.18  ms)  and  tenth  (0.72  ms)  time  channels.  The  plots  were  obtained  by 
stitching  together  the  two  data  sets.  This  was  achieved  by  averaging  the  data  collected  on  duplicate  stations.  For 
presentation  pinposes,  the  data  in  each  contour  plot  were  thresholded  from  below  at  a  1  mV  level  and  truncated 
above  at  levels  indieated  in  the  plot  titles. 
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Figure  1 1 :  4  m  x  4  m  area  on  which  a  series  of  EM63  surveys  were  taken  over  different  targets  at  several  depths  and  orientations. 


(a)  Photograph  of  the  wooden  target  holder  for  EM63  mea-  (b)  Measurement  procedure  when  using  the  jig. 

surements.  A  lOSmm  projectile  is  standing  beside  the  jig. 


!  Figure  12:  Setup  for  EM63  decay  measurements 

I 

I 

1 
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Voltag*  (mv) 


Figure  13:  TEM  soundings  measured  using  jig.  Unfilled  symbols  (e.g.  ’o’)  indicate  negative  data  that  has  its  absolute  value 
plotted.  Predicted  response  is  obtained  by  fitting  the  measurement  with  Equation  24. 
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(a)  Photograph  of  setup  used  to  measure  response  of  smaller  scrap  targets. 


(b)  Voltage  decay  measurements  for  scrap  targets  1  to  4.  (c)  Voltage  decay  measurements  for  scrap  targets  5  to  8 

Figure  14:  Setup  and  sample  plotting  response  measurements  of  small  non-UXO  items.  Unfilled  symbols  (e.g.  ’o’)  indicate 
negative  data  that  has  its  absolute  value  plotted. 
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(a) 


TlmeChamell:  0.18  ntMc 


(b) 


TImaChannail:  0.18fn«ec 


Figure  15:  Measured  response  of  scrap  targets  1  to  8  along  a  line.  Panels  (a)  and  (b)  plot  the  response  along  a  line  for  the  first 
time  channel  {t  =  0.18  msec).  Panels  (c)  and  (d)  plot  the  response  along  a  line  for  the  tenth  time  channel  (t  =  0.72  msec). 
Each  target  is  located  at  station  0  cm. 
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(channel  10).  Truncate  =100 


2  Verification  of  the  Dipole  Model 


In  Part  I  of  this  report  an  approximate  forward  model  was  proposed  for  the  TEM  response  of  a  buried  axi-symmetric 
metallic  target.  In  this  forward  model,  the  secondary  field  is  approximated  by  a  pair  of  orthogonal  and  indepen¬ 
dently  decaying  dipoles,  whose  strengths  are  proportional  to  the  projection  of  the  primary  field  onto  their  direction. 
The  decay  of  each  dipole  moment  is  governed  by  the  magnetic  polarization  tensor  M.  The  magnetic  polarizability 
tensor  is  independent  of  transmitter/receiver/target  geometry  and  is  a  function  of  the  physical  characteristics  of  the 
target  alone.  In  our  previous  work  we  outlined  a  technique  for  exploiting  M  as  a  tool  for  characterizing  the  shape 
of  target  as  either  rod-like  or  plate-fike,  and  if  the  target  is  ferrous  or  nonferrous  (Pasion  1999).  These  results  can 
only  be  applied  with  confidence  for  measurement  configurations  where  the  forward  modeUmg  is  applicable. 


This  section  contains  a  series  of  tests  designed  to  verify  that  the  magnetic  polarization  tensor  is  indeed  independent 
of  transmitter/receiver/taiget  geometry.  Each  test  follows  the  same  procedure.  First,  the  target  is  measured  in  two 
orientations:  with  the  axis  of  symmetry  parallel  to  the  primary  field  and  then  perpendicular  to  the  primary  field. 
These  measurements  allow  us  to  extract  the  target’s  two  characteristic  decay  curves  L  i  [t)  and  L2  [t)  that  define 
the  magnetic  polarization  tensor.  Equipped  with  we  can  then  predict  the  TEM  response  for  various  locations 
and  orientations  and  compare  these  results  with  measured  data.  In  this  section  we  evaluate  the  accuracy  of  these 
predictions. 


Obtaining  the  Polarizability  Tensor 

Constructing  the  polarizability  tensor  requires  obtaining  the  dipole  decay  functions  Li  (t)  and  L2  {t).  The 
decay  functions  can  be  isolated  by  making  two  measurements:  (1)  primary  field  B  ^  parallel  to  the  z'  axis  of 
symmetry,  and  (2)  primary  field  perpendicular  to  the  z'  axis  of  symmetry. 


Figure  17  illustrates  the  arrangement  of  the  EM63  and  target  that  we  used  to  obtain  the  decay  functions.  In  this 
geometry,  the  measured  voltage  is  then: 


V  {t)  =  K 


2  [Z/i  (t)  cos^  0  +  L2  (t)  sin^  0] 


(24) 


where  Z  is  the  distance  between  the  center  of  the  receiver  loop  and  the  center  of  the  target,  and  k  is  a  constant 
that  depends  on  the  size  of  the  receiver  and  transmitter  loops,  the  munber  of  turns  in  each  loop,  and  the  transmitter 
current. 


When  the  target’s  z'  axis  of  symmetry  is  parallel  to  the  inducing  field  (9  =  0  degrees.  Figure  17a),  only  the  m  1  (t) 
dipole  is  excited,  and  the  measured  voltage  is 


T/ll  it)  = 


B^iZ) 


2Li  (t)  = 


2k, 


B^iZ) 


A;i(t-bai)  e 


(25) 


and  when  the  target’s  z'  axis  of  symmetry  is  perpendicular  to  the  inducing  field  (6  =  90  deg.  Figure  17b),  then 
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only  the  m2  (t)  dipole  is  excited  and  the  observed  voltage  is: 


(t)  =  K 


B^jZ) 


2L2  (t)  = 


2k 


BP{Z) 


k-z  (t  +  02) 


-02  g-th2 


(26) 


The  voltage  curves  recorded  by  the  EM63  for  the  parallel  and  perpendicular  responses  can  be  obtained  by  using 
a  scaled  down  version  of  the  inversion  algorithm  to  recover  the  decay  parameters  (Pasion  1999).  Figures  17c 
and  d  give  the  result  of  this  procedure  applied  to  a  105  mm  projectile  and  scrap  14.  As  expected  for  a  rod-like 
target  the  Li  component  of  the  polarization  tensor  has  a  greater  magnitude  than  the  L  2  component.  Figure  18  has 
the  recovered  Li  (t)  and  £2  (t)  for  the  targets  placed  in  the  jig.  Several  of  the  curves  have  been  extended  by  a 
“dash-dot”  line  that  indicates  an  extrapolated  portion  of  L 1  (t)  and  £2  (^)- 

The  accuracy  of  this  procedure  to  obtain  £  1  (t)  and  £2  (t)  will  depend  on  experimental  error  and  model  error. 
Possible  experimental  errors  include:  (1)  inaccurate  measurement  of  Z;  (2)  inaccurate  placement  of  the  target 
beneath  the  receiver  loop,  i.e.,  the  center  of  the  target  must  be  placed  along  the  vertical  axis  passing  through  the 
center  of  the  loop;  (3)  tilting  of  the  EM63  trailer  such  that  the  primary  field  isn’t  vertical.  Modelling  errors  describe 
instances  where  assumptions  of  the  forward  model  are  violated,  such  as:  (1)  uniformity  of  the  primary  field  in  the 
volume  of  the  target;  (2)  representing  the  response  as  a  point  dipole;  and  (3)  absence  of  fore-aft  symmetry. 


Test  1 :  Reproducing  the  Time  Sounding  at  Intermediate  Orientation 

Once  we  determine  the  parameters  of  the  magnetic  polarization  tensor,  we  can  forward  model  the  parameters  to 
obtain  the  TEM  response  for  any  location,  depth,  and  orientation.  In  this  test  we  investigate  how  accurately  we 
could  predict  the  decay  of  the  secondary  field  at  an  arbitrary  orientation  for  a  recovered  magnetic  polarization 
tensor. 

Equation  24  describes  the  measured  voltage  in  a  receiver  coil  directly  above  a  target  illuminated  by  a  purely 
vertical  primary  field.  Using  Equation  24  and  the  decay  functions  £  1  (t)  and  £2  (t)  obtained  in  the  analysis  of 
the  previous  section,  we  predict  the  voltage  response  and  compare  it  to  measurements  using  the  EM63.  The  set 
of  measurements  we  use  for  comparison  were  those  obtained  using  the  target  holder.  Figures  19  and  20  compare 
the  measured  responses  at  different  angles  and  the  response  predicted  by  the  forward  model  for  an  81  mm  mortar 
without  fins.  Figure  19  has  the  measured  voltage  curves  at  different  heights  from  the  sensor.  £  1  (t)  and  £2  (t), 
obtained  by  fitting  these  curves  with  Equations  25  and  26,  are  then  forward  modelled  and  plotted  in  Figure  20.  The 
top  two  panels  of  Figure  20  demonstrate  the  procedure  on  an  81  mm  mortar  (without  fins)  located  approximately 
55  cm  beneath  the  receiver  loop.  At  this  distance  (equivalent  to  approximately  10  cm  below  the  surface)  we  see 
that  the  representation  is  only  moderately  good  at  reproducing  the  data  at  the  different  angles.  When  we  repeat  the 
procedure  for  data  collected  approximately  100cm  beneath  the  receiver,  we  see  that  the  model  does  abetter  job  of 
predicting  the  data.  This  is  not  surprising,  since  we  would  expect  the  modeUing  of  a  compact  metallic  object  as  a 
dipole  to  become  more  applicable  as  we  move  further  fi'om  the  source/receiver  loop.  Figure  21  demonstrates  the 
procedure  for  predicting  the  measured  voltage  decay  curves  for  Scrap  5  using  the  dipole  model. 


Test  2:  Reprociucing  the  Spatial  Behavior 

In  our  second  test,  we  focus  on  how  accurately  we  could  predict  the  spatial  response  of  a  target  using  the  dipole 
model.  For  this  investigation  we  first  use  a  target’s  magnetic  polarization  tensor,  obtained  in  the  manner  described 


Chapter  2  Verification  of  the  Dipole  Model 


39 


in  the  previous  section,  to  predict  the  secondary  field  over  a  survey  line  that  passes  directly  over  the  target.  In  each 
case,  the  survey  line  is  co-aligned  with  the  target.  This  predicted  response  is  then  compared  with  the  measured 
response. 


In  Figures  22  to  30,  the  predicted  and  measured  responses  along  a  line  are  plotted  for  several  targets  in  three 
orientations;  vertical,  horizontal,  and  at  an  intermediate  dip  angle.  At  intermediate  angles  the  target  dips  toward 
the  end  of  the  line  (i.e.,  dips  downward  to  the  right  in  Figures  22  to  30).  Figures  22  to  30  indicate  the  model  is 
successful  in  predicting  the  response  along  each  survey  line. 
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0.180  1.0  10.0  25.14 

Time  (ms) 


Figure  19:  Measured  voltage  curves  for  a  81  mm  Mortar  without  6ns  at  different  heights  from  the  sensor.  Un611ed  symbols  (e.g. 

indicate  negadve  data  that  have  absolute  value  plotted.  These  curves  were  inverted  to  obtain  Li  (f)  and  L2  (t).  Clearly, 
the  decay  nature  of  the  signal  changes  when  the  mortar  is  brought  close  to  the  sensor.  The  resulhng  Li  (t)  and  L2  (t)  curves 
are  used  to  predict  the  decay  at  intermediate  angles  in  Figure  20. 


Chapter  2  Verification  of  the  Dipole  Model 


43 


Z  ~  55  cm 


Z  ~  55  cm 


Time  (ms) 


Time  (ms) 


Z- 100  cm  Z- 100  cm 


0.180  1.0  10.0  25.14  0.180  1.0  10.0  25.14 

Time  (ms)  Time  (ms) 


Figure  20:  Predicting  the  measured  voltage  decay  curves  for  an  81  mm  mortar  without  fins  at  different  orientations  9.  Unfilled 
symbols  (e.g.,  “□’)  indicate  negative  data  that  have  absolute  value  plotted.  The  L\  (<)  and  L2  {t)  decay  curves  recovered  from 
measurements  with  Z=60  cm  and  Z=52  cm  (Figure  19),  respectively,  are  used  to  predict  the  voltage  curves  in  the  top  two 
panels.  The  Li  (f)  and  L2  (f)  decay  curves  recovered  from  measurements  with  Z=102  cm  and  Z=95  cm,  respectively,  are 
used  to  predict  the  voltage  curves  in  the  bottom  two  panels.  The  prediction  of  the  measured  curves  is  more  successful  when  the 
target  is  farther  away  from  the  sensors. 
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Figure  24:  60  mm  mortar:  0  =  40  deg,  Z=15  cm. 
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Figure  29:  Stokes  mortar:  0  =90  deg  (horizontal),  Z=73  cm. 
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Figure  30;  Stokes  mortar;  0  =  54  deg,  Z=73  cm. 
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3  Inversion  of  TEM  Data  Collected  at  the  ERDC 
Test  Site 


In  Part  I  of  this  report  we  described  a  method  of  inverting  TEM  data  for  the  13  parameters  of  the  Two-Dipole 
model.  We  formulated  a  nonlinear  least-squares  problem  that  involved  minimizing 

^H  =  ^||Wd(F[m]-d'>*>*)|P  (27) 

where  F[m]  is  the  forward  modelled  data,  d  is  the  observed  data,  and  $  is  the  least  squares  objective  function 
that  measures  how  closely  our  predicted  data  match  the  observed  data.  W  ^  is  the  data  weighting  matrix 


1 

pdi  -t-  e 


where  p  is  typically  a  percentage  and  e  is  a  constant  that  characterizes  ambient  noise. 


(28) 


In  Part  I  we  performed  two  transformations  to  the  model:  (1)  a  simple  linear  scaling  to  the  model  m  to  ensure 
that  each  component  of  m  is  approximately  unity,  and  (2)  replacing  each  of  the  time  decay  parameters  m  ,  by  a 
squared  variable  Wi  such  that  m*  =  x?  to  ensure  positivity  of  the  decay  parameters.  For  the  inversions  in  Part  n 
of  this  report  we  retain  the  Unear  scaling  of  the  model,  but  we  no  longer  solve  the  square- variable  unconstrained 
problem.  Instead  we  reformulate  the  inversion  so  that  it  can  be  solved  by  a  constrained  optimization  procedure.  We 
use  a  Projected  BFGS  algorithm,  described  in  “Iterative  Methods  for  Optimization”  (Kelley  1999),  to  minimize 
the  objective  function.  In  this  algorithm,  we  are  required  to  supply  a  starting  model  as  well  as  upper  and  lower 
constraints  for  each  variable. 


Selection  of  Starting  Parameters 

The  first  step  of  this  inversion  procedure  is  to  make  a  starting  guess  for  the  13  model  parameters: 

m  =  [X,Y,Z,(l),e,ki,ai,Pi,'yi,k2,a2,P2,72]  (29) 

The  success  of  a  local  inversion  procedure,  as  well  as  the  rate  of  convergence  toward  a  solution  is  dependent  on  the 
quality  of  the  initial  guess  (Bard  1974).  The  following  section  will  discuss  a  number  of  simple  data  preprocessing 
strategies  that  enable  us  to  make  a  reasonable  initial  guess  of  the  above  parameters. 

Location  on  Survey  {X,  Y) 

A  natural  idea  for  determining  the  location  for  a  target  would  be  to  find  the  location  at  which  the  maximum  signal 
occurs.  This  technique  would  produce  the  exact  location  for  a  spherical  target,  where  the  plot  of  the  vertical 
(z)  secondary  field  produces  a  circular  “buUs-eye”  anomaly  with  a  maximum  located  directly  over  the  target 
(Figure  31a).  However  rod-fike  and  plate-like  targets  can  produce  non-synunetric  anomalies  with  multiple  local 
maxima  (Figure  31b).  In  these  cases,  search  for  the  maximum  signal  wiU  give  an  inaccurate  estimation  of  target 
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location.  To  account  for  these  possible  anomaly  types,  we  use  an  approximation  of  the  first  spatial  moment  of  the 
data  to  estimate  the  location  (A'^  ,  y  «**)  on  the  survey: 


Xe^t 


E^l  Vj  ih)Xi 


and 


est  _  Vj  jtl)  Vi 


Y 


EliViih) 


(30) 


where  there  are  N  total  stations  in  the  survey  and  Fj  (fi)  is  the  voltage  measured  at  the  F'*  station  at  location 
{xi,yi),  and  at  the  first  time  channel.  Figure  31  compares  this  technique  to  simply  using  the  maximum  value. 


Depth  of  target  (Z) 

The  full  width  at  half  maximum  (FWHM)  is  a  standard  technique  for  estimating  the  depth  of  a  target.  In  general, 
the  plotted  anomaly  of  nonspherical  targets  wiU  not  be  circularly  symmetric.  In  these  cases,  it  is  not  clear  which 
profile  to  use  for  calculating  the  FWHM.  Therefore  we  modify  the  FWHM  technique  by  plotting  points  of  the 
secondary  field  anomaly  at  an  early  time  channel  that  are  greater  than  one-half  of  the  maYimnm  signal.  The 
estimate  for  the  depth  of  the  target  is  then: 


Estimated  Depth  = 


(31) 


where  A  is  the  area  of  the  portion  of  the  secondary  field  anomaly  greater  than  one-half  of  the  maximiun  signal. 


Orientation  {(p,  6) 

The  orientation  angles  (}>  and  6  are  not  estimated  prior  to  the  inversion.  The  starting  values  for  </>  and  9  set  to  45 
deg. 


Decay  Parameters  {h,  /5i,  7^) 

The  decay  parameters  contained  within  the  decay  functions  L  i  (t)  and  L2  (t)  are  estimated  in  the  same  manner 
described  by  Pasion  (1999).  The  decay  parameters  are  estimated  by  taking  the  estimate  of  the  target  location  and 
extracting  the  decay  curve  at  that  location.  The  decay  curve  is  then  fit  with  the  function: 

at)  =  K^-^-T-2k{a  +  tfe*/^  (32) 

where  is  the  estimated  depth.  The  starting  decay  parameters  for  both  dipoles  1  and  2  are  then  ( k,  a,  0, 7). 
Equation  32  is  obtained  by  setting  L 1  —  L2  in  Equation  24. 


Selection  of  Upper  ar\6  Lower  Limits 

The  nonlinear  least-squares  objective  function  (Equation  27)  is  minimized  using  an  optimization  algorithm  that 
incorporates  box  constraints.  Box  constraints  are  the  lower  and  upper  limits  that  we  place  on  the  model  parameters 
that  represent  the  minimum  and  maximum  possible  values  of  the  model  parameters.  The  upper  and  lower  limits 
used  for  the  model  parameters  are  summarized  in  Table  3. 


Location  (A",  Y)  and  Depth  (Z) 

For  the  inversion  results  presented  here,  conservative  constraints  were  apphed.  The  location  was  constrained  to  a 
1  m  square  box  centered  at  (A",  Y)  =  (2  m,2  m),  and  the  depth  was  constrained  between  Z  low  —  0.43  m,  which  is 
the  level  of  the  ground,  and  Z^igh  =  1.5  m. 

Chapter  3  Inversion  of  TEM  Data  Collected  at  the  ERDC  Test  Site 


54 


Orientation  (0, 9) 

Since  we  do  not  estimate  the  orientation  of  the  target  based  on  the  plotted  anomaly,  the  upper  and  lower  limits  of  cj) 
and  6  are  determined  to  be  the  limits  that  allow  for  all  possible  orientations  of  the  target.  Therefore,  -90  <  ^  <  90 

and  0  <  <  180. 


Decay  Parameters  (/ci,  a,,  A,  7*) 

In  verifying  the  applicability  of  the  approximate  forward  model  (Chapter  2,  Part  H),  the  decay  parameters  were 
recovered  for  different  scrap  and  UXO  targets.  From  this  analysis,  and  the  analysis  reported  by  Pasion  (1999), 
upper  and  lower  limits  for  the  decay  parameters  were  chosen.  The  limits  are  listed  in  Table  3b. 


Application  to  Field  Data  Sets 

In  this  text  we  present  the  results  of  applying  the  inversion  algorithm  to  five  targets:  155  mm,  105  mm,  81  mm,  60 
tntn,  and  Stokes  mortar.  A  list  of  inversion  results  presented  in  this  report  can  be  foimd  in  Table  4,  where  the  target 
type,  depth,  and  orientation  for  each  data  set  inverted  in  this  section  are  given.  Each  inverted  data  set  consisted  of 
soundings  collected  on  a  2  x  2  m  grid  centered  on  the  target,  containing  five  lines  running  North-South  separated 
at  50  cm  line  spacing,  with  stations  located  at  5  cm  intervals  along  each  line.  Inversion  results  can  be  found  in 
Figures  32  to  48.  In  Panel  a  of  each  figure  is  a  plan  view  comparison  of  the  observed  and  predicted  data  for  time 
channels  1,  8, 15,  and  23.  Panels  b  compare  the  predicted  and  observed  decay  at  three  stations:  Line  1  m.  Station 
1  m;  Line  2  m.  Station  2  m;  and  Line  2.5  m.  Station  2.5  m.  The  recovered  location  and  orientation  are  listed  in 
Panels  c,  and  the  recovered  decay  parameters  and  diagnostics  applied  to  those  parameters  are  reported  in  Panels  d. 
The  algorithm  was  successful  in  identifying  each  target  as  tod-like  and  permeable  for  aU  the  cases  presented  here. 
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Figure  31:  Finding  the  location  of  a  (a)  sphere  and  (b)  horizontal  stokes  mortar.  The  location  obtained  by  using  the  maximum 
of  the  signal  is  indicated  by  the  star,  and  the  location  obtained  by  using  the  “center-of  mass”  is  indicated  by  the  triangle.  The 
center-of-mass  is  much  better  at  determining  the  location  of  the  target  when  the  plotted  anomaly  does  not  have  the  appearance 
of  a  “bullseye”  (e.g.,  a  sphere).  Panels  c  and  d  plot  those  data  points  that  are  greater  than  half  the  maximum  signal.  Panel  d 
indicate  that  the  traditional  “FuU-Width-at-Half-Maximum”  would  be  difficult  to  utilize  for  non-“bullseye”  anomalies. 


Table  3:  Upper  and  lower  limits  used  in  inversion  examples  in  Part  n  of  this  report 


m-i 

m7 

X(m) 

1.5 

2.5 

Y(m) 

1.5 

2.5 

Z  (m) 

0.43 

1.5 

(j)  (degrees) 

-90 

90 

9  (degrees) 

0 

180 

rrii 

up 

m7 

ki 

0.1 

200 

a* 

0.001 

0.02 

Pi 

0.3 

1.5 

7* 

2 

30 

a.  Location  and  Orientation 


b.  Decay  Parameters 
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Table  4:  Targets  used  in  field  set  examples 


Target 

Z  (cm) 

6  (degrees) 

Figure  No. 

155  mm 

98 

0 

32 

155  mm 

97 

52 

33 

155  mm 

91 

90 

34 

105  mm 

86 

0 

35 

105  mm 

87 

90 

36 

105  mm 

89 

59 

37 

105  mm 

89 

40 

38 

Stokes  Mortar 

73 

54 

39 

Stokes  Mortar 

73 

0 

40 

Stokes  Mortar 

73 

90 

41 

81  mm 

73 

0 

42 

81  mm 

74 

43 

43 

81  mm 

Ti 

90 

44 

60  mm 

58 

57 

45 

60  mm 

58 

90 

46 

60  mm 

59 

0  (nose  up) 

47 

60  mm 

59 

0  (nose  down) 

48 
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Figure  34:  Recovered  parameters  for  the  field  data  inversion  of  a  155  mm  projectile. 


105  mm:  0  =  0  degrees  (vertical),  Z  =  86  cm 
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105  mm:  (9  =  90  degrees  (horizontal), 
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Figure  37:  Recovered  parameters  for  the  field  data  inversion  of  a  105  mm  projectile. 


Observed  Predicted 


Figure  38:  Recovered  parameters  for  the  field  data  inversion  of  a  105  mm  projectile. 


Figure  39:  Recovered  parameters  for  the  field  data  inversion  of  a  stokes  mortar. 


Stokes  Mortar;  9  =  0  degrees  (vertical),  Z  =73  cm 


Figure  40;  Recovered  parameters  for  the  field  data  inversion  of  a  stokes  mortar. 


Figure  41:  Recovered  parameters  for  the  field  data  inversion  of  a  stokes  mortar. 


81  mm:  0  =  0  degrees  (vertical). 


Figure  42:  Recovered  parameters  for  the  field  data  inversion  of  an  81  mm  mortar  without  fins. 


Observed  Predicted 


Figure  43:  Recovered  parameters  for  the  field  data  inversion  of  an  81  mm  mortar  without  fins. 


81  mm:  9  =  90  degrees  (horizontal),  Z  =  73  cm 
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Figure  44:  Recovered  parameters  for  the  field  data  inversion  of  an  81  mm  mortar  without  fins. 


Figure  45:  Recovered  parameters  for  the  field  data  inversion  of  a  60  mm  mortar. 


mm  mortar. 


Observed  Predicted 


60  imn:  0  =  0  degrees  (vertical,  nose  down),  Z  =  59  cm 
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Figure  48:  Recovered  parameters  for  the  field  data  inversion  of  a  60  mm  mortar. 
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A  Target  Descriptions 


Target 

Mass 

(kg) 

Length 

(cm) 

Width 

/Thickness  (cm) 

155  mm 

60 

15.5 

105  mm 

14.05 

21 

10.5 

Stokes  Mortar 

3.55 

36 

16 

81  mm  (no  fins) 

3.3 

26 

8.1 

Rusted  Mortar 

2.85 

36 

29 

60  mm 

1.352 

26 

6 

37  mm 

0.839 

11.4 

3.6 

Disk 

0.88 

8  (diameter) 

2.2  thick 

Scrap  1 

0.08 

13  (diameter) 

~2  thick 

Scrap  2 

0.026 

6.4 

6.3 

Scrap  3 

0.055 

11 

6.4 

Scrap  4 

0.091 

19.4 

2.5 

Scrap  5 

0.172 

16.5 

3 

Scrap  6 

0.256 

21 

4.8 

Scrap  7 

0.069 

12.5 

2 

Scrap  8 

0.075 

7 

3 

Scrap  9 

19 

9 

Scrap  10 

0.032 

5 

3.5 

Scrap  12 

0.186 

9.6 

4.4 

Rocket  Fins 
(Scrap  13) 

0.83 

21.4 

6.5 

Blown  Mortar 
(Scrap  14) 

0.939 

16 

7 

Scrap  15 

0.431 

16 

5.6 

A  Target  Descriptions 


Rusted 

Mortar 


Stokes 

Mortar 


155imn 


Figure  49:  UXO  targets  analysed. 


B  Recovered  Polarization  Dyadic  for  Various 
Targets 


The  magnetic  polarization  dyadic  for  the  TEM  problem  is  defined  as 


M  = 


'L2  [t) 
0 
0 


0 

L2  {t) 
0 


0 

0 

Li  (t) 


where 


Li  {t)  =  h  + 

L2  (t)  =  ^2  (t  +  02) 


(33) 


(34) 


The  recovered  magnetic  polarization  dyadic  acquired  for  a  number  of  targets  are  reported  here.  The  procedure 
followed  to  obtained  these  results  is  described  in  Chapter  2  of  Part  n  of  this  report.  A  photo  and  description  of 
each  target  can  be  found  in  Appendix  A. 


(b)  Decay  Parameters 


mi 

Dipole  1 

Dipole  2 

k 

69.6 

20.1 

a 

0.002 

0.005 

p 

0.64 

1.08 

7 

20.4 

7.59 

Figure  51:  105  mm  Projectile 
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Rusted  Mortar 


(b)  Decay  Parameters 


Scrap  14 


(b)  Decay  Parameters 


rrii 

Dipole  1 

Dipole  2 

k 

1.78 

0.89 

a 

0.02 

0.02 

p 

0.80 

1.11 

7 

8.99 

2.14 

Figure  58:  Scrap  14 


(b)  Decay  Parameters 


k 

1.26 

na 

a 

na 

p 

na 

7 

2.08 

na 

Figure  59:  Scrap  5.  The  response  of  the  m2  dipole  did  not  generate  a  response  above  the  noise  level  of  the  EM63. 


(b)  Decay  Parameters 


rrii 

Dipole  1 

Dipole  2 

k 

0.23 

1.45 

a 

0.001 

0.02 

p 

1.44 

0.99 

7 

2.57 

2.00 

Figure  60:  Steel  disk 
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